#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "common.h"

//
//  benchmarking program
//
int main( int argc, char **argv )
{    
	if( find_option( argc, argv, "-h" ) >= 0 )
	{
		printf( "Options:\n" );
		printf( "-h to see this help\n" );
		printf( "-n <int> to set the number of particles\n" );
		printf( "-o <filename> to specify the output file name\n" );
		return 0;
	}

	int n = read_int( argc, argv, "-n", 1000 );

	char *savename = read_string( argc, argv, "-o", NULL );

	FILE *fsave = savename ? fopen( savename, "w" ) : NULL;
	particle_t *particles = (particle_t*) malloc( n * sizeof(particle_t) );
	set_size( n );
	init_particles( n, particles );

	//-----------------------------------------------------

	// Creates the main box, this is box number 0, the root
	box world;
	init_world(&world);
	
	// Adds the particles to the world box
	for(int i = 0; i < n; i++ )
		world.add_new(i);
	
	// Calculates what the depth level should be
	// The depth level is the number of layers our simulation will have
	int depth_level = ceil(log10(n));
	printf("depth_level = %d \n", depth_level);

	// Each layer is a 2 dimensional array containing boxes
	// layer 0 has 1 box
	// layer 1 has 4 boxes
	// layer 2 has 8 boxes
	// etc...
	box ***boxes;	// This will be the array containing all of the boxes/simulation domains
	boxes = new box**[depth_level + 1];	// Creates the level 0 multi dimensional array

	// box_list is a list each thread has, in the case of serial there is only 1 box_list
	// box_list contains the list of all the boxes, on each layer, that fmm must go through and calculate
	// This is done here so that fmm can just do calculations instead of figuring out in which
	// boxes to do the calculation
	vector < vector < pair < int, int > > > box_list;
	for(int i = 0; i <= depth_level; i++)
	{
		vector < pair<int, int > > dummy_v;
		box_list.push_back(dummy_v);
	}
	
	// Initializes the 2 dimensional arrays, creates the boxes
	initializations(boxes, depth_level);
	// Links the root boxes 4 children to it and sends them their particles
	split(world, particles, boxes, 0, 0);
	// Links the 4 level 1 children with their own children
	// and finishes linking all the other boxes with their children
	// also, finishes distributing all the particles to each box accordingly to their coordinates
	build_tree(boxes, particles, 1, depth_level, 0, 0, 1, 1);

	// All we need to do now is to 

	// This part wanders through the level 2 boxes this process got assigned to itself
	// It calls the function find_bottom_box(), this function goes to the bottom layer
	// children of this box and add them to this process's bottom_boxes list.
	int tcx = 0;
	int tcy = 0;
	// 16 is used because level 2 depth has 16 boxes
	// this is a good depth to use because it allows us to have 16 threads
	// from here if it is pthreads, openmp or mpi and it is not always necessary to 
	// go to level 3 on smaller simulations
	for(int nn = 0; nn < 16; nn++)
	{
		if(tcy > 3)
		{
			tcy = 0;
			tcx++;
		}
		find_bottom_box(boxes[2][tcx][tcy], boxes, box_list, depth_level, tcx, tcy);
		tcy++;
	}

	
	//
	//  simulate a number of time steps
	//  
	int i;
	double simulation_time = read_timer( );
	for( int step = 0; step < NSTEPS; step++ )
	{
		// Clears all of the forces applied on the particles
		for( int i = 0; i < n; i++ )
			particles[i].ax = particles[i].ay = 0;
		
		// Calls Fast Multipole Method
		fmm(boxes, particles, depth_level, box_list);
		
		// Moves the particles
		for( int i = 0; i < n; i++ ) 
			move( particles[i] );
		
		// The particles have been moved, now the get resorted to each box
		sort_elements(world, particles, boxes, 0, 0);
		assign_elements(boxes, particles, 1, depth_level, 0, 0, 1, 1);
		
		// save if necessary
		if( fsave && (step%SAVEFREQ) == 0 )
			save( fsave, n, particles );
	}
	
	simulation_time = read_timer( ) - simulation_time;

	printf( "n = %d, simulation time = %g seconds\n", n, simulation_time );

	free( particles );
	if( fsave )
	fclose( fsave );

	return 0;
}
